"Glassy" Relaxation in Catalytic Reaction Networks 
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Relaxation dynamics in reversible catalytic reaction networks is studied, revealing two salient 
behaviors that are reminiscent of glassy behavior: slow relaxation with log(time) dependence of the 
correlation function, and emergence of a few plateaus in the relaxation. The former is explained 
by the eigenvalue distribution of a Jacobian matrix around the equilibrium state that follows the 
distribution of kinetic coefficients of reactions. The latter is associated with kinetic constraints, 
rather than metastable states, and is due to the deficiency of catalysts for chemicals in excess 
and negative correlation between the two chemical species. Examples are given, and generality is 
discussed. 

PACS numbers: 87.16.Yc, 82.39.Rt, 05.40.-a 



Cells are usually not in thermal equilibrium, and bi- 
ological functions are believed to operate under non- 
equilibrium conditions. The relevance of non-equilibrium 
conditions to pattern formation has been discussed for 
decades 0] since the pioneering work of Schr6dinger[H. 
In contrast to physics and chemistry, however, such non- 
equilibrium conditions are not imposed externally but 
have to be sustained by a biological system itself. This 
sustainment might then suggest the existence of some 
bootstrapping process in which biochemical reactions un- 
der non-equilibrium conditions could suppress relaxation 
to equilibrium. Even though this argument may be too 
naive for currently known living organisms that adopt 
advanced mechanisms using cell membranes, it is never- 
theless important in considering the origin of life. 

In physics, the reluctance to relax to equilibrium has 
been studied in glass, and a certain coniplex free energy 
landscape structure has been elucidated [1, 0,11, Hi- As an 
alternative to such structural studies, kinetic mechanisms 
to suppress the relaxation have recently been proposed 
0, 'Kinetically constrained models' have gath- 

ered much attention 11 , 13, l3| , where the relaxation to 
equilibrium is slowed down due to a kinetic bottleneck. 
In the present Letter, we show that, in a system with a 
catalytic reaction network, relaxation to thermal equilib- 
rium is generally slowed down due to a kinetic constraint. 

We consider a network of reactions consisting of M 
chemical components {Xi, z = 1, • ■ • , M), each of which is 
catalyzed by one of the M components. Transformation 
between chemicals Xi and Xj is catalyzed by Xc{i,j), 



X, + X 



(1) 



sume that all chemical species are percolated to any other 
through these reactions. The system is closed, without 
inflow of chemicals or energy from the outside. Note 
that the number of molecules, accordingly Xi = S, is 
conserved by the above reactions, where Xi is the concen- 
tration of each chemical species i. 

To assure the relaxation to thermal equilibrium, the 
ratio of forward to backward reactions is set so that it 
satisfies the detailed balance condition. It is satisfied 
by allocating energy Ei to each molecular species, and 
setting the ratio of forward (fcij) to backward (fcj.i) re- 
actions in eq. (1) to ki,j/kj^i = exp(— /3(£'j — Ei)), where 
(3 is the inverse temperature. As a result, the equilib- 
rium concentration xi'' satisfies xi'' = sexp(—PEi) with 

s = s{j:,eM-m))-'- 

Here we take a continuum description, so that the dy- 
namics of the concentration is given by the rate equation 



(2) 



The reaction network consists of the above reactions, 
with the total number of reactions G > M. We as- 



with ki_j = min{l, exp(— /3(£^j — ii^j))}, and Con{i, j] c) — 
Con{j, i;c) = 1 if there is a reaction path, as in eq. (1), 
and otherwise [l^. Note that eq. (2) has a unique stable 
fixed point attractor xi'', without any metastable states. 
We assume that the energy Ei is distributed uniformly, as 
jj£ (e is a constant) [lH|. The network C'on{i,j; c) is cho- 
sen randomly by setting the average number of paths for 
each chemical K — 2G/M. As an example of a typical re- 
laxation course, we set an initial concentration with equal 
distribution over all chemicals, i.e., the high-temperature 
limit (corresponding to = 0), and study the evolution 
under given (3. 

In Fig. 1, we give examples of the relaxation time 
course for four sets of networks (M ~ 24, K — 8), 
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FIG. 1: (a)(b) Relaxation time course for four sets of networks 
(M = 24, K = 8) for several /3. 
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FIG. 2: Relaxation time as a function of /3 for the sample 
reaction networks in Fig. l(a)(b). 



where we plot the deviation from equilibrium concen- 



tration defined by C{t) = {J2i(^i(^) 



)}/{^-(xi(0) — x^"^)^}. We note two salient behaviors 
when (3 is sufficiently larger than Pc ~ 3/e, which is the 
inverse of the average difference between energy levels. 
First, there exists overall log(t) relaxation, in contrast to 
exponential relaxation for small /3. Second, there are sev- 
eral plateaus in the relaxation course. The logarithmic 
relaxation at large /3 is generally observed, independently 
of the networks or K. Existence of plateaus is also uni- 
versal. The number of plateaus depends on each network 
(Fig. 1(b)); generally, the number decreases as K in- 
creases. 

The integrated relaxation time r =< \C{t)\dt > is 
plotted as a function of /3 for several M and K in Fig. 
2, where < ■ • • > is the average over networks with given 
M and K. For the small /? regime in which C{t) decays 
exponentially, t follows exp(/3/ /3c). It is the inverse of the 
average reaction rate to increase the energy, which gives 
the order of the relaxation time[l3l. For large giving 
log(t) relaxation, r follows exp{Rf3e) with R approaching 
a larger constant with the increase in (3. R increases with 
the increase (decrease) in M (K), respectively. 



The log(t) relaxation with plateaus is often observed in 
glass theory and experiments (HI. [lH . In the present case, 
these relaxation characteristics are partially explained by 
a rough estimate of the eigenvalue distribution in linear 
stability analysis. Consider deviation from the equilib- 
rium concentration as Xi(t) — x^^ + Sxi{t), where the 
equilibrium concentration x^' — s exp{—PEi) is the fixed 
point solution of eq. (2). By linearizing with Sxi{t) {i = 
1, • • • , M), we get (5x(t) = J5x(t) with the Jacobi matrix 
J computed straightforwardly. For large /3, Jij for i > j, 
given by Con{i,j; c')x'^Je~^^^'^~^^\ is much smaller than 
that for i < j, Con{i,j;c")x'^?,. If the former terms are 
neglected, the above J is a triangular matrix, so that the 
eigenvalues Xi of J are given by diagonal elements Ji i = 
- Ej<. Con[i,r, c")4? -E,>, Con{^,J■ d)xlU^^^^-^^\ 
whose distribution has similar dependence to that of 
exp(— /JiJfe), for large (i. This is also true for the neglected 
off-diagonal terms. Hence, it is expected that the distri- 
bution of the eigenvalues \i is similar to the distribution 
of cxp{—(3Ek), for large f3 (besides the null eigenvalue 
Ao = corresponding to the equilibrium distribution). 
In fact, numerical diagonalization of the Jacobian ma- 
trix supports this estimate of eigenvalue distribution. By 
using this linear approximation and the correspondence 
of the eigenvalue with exp(— /3i?), C(t) is approximated 
by D{E)a{E)exp{-e-'^'^t)dE, with the distribution 
of energy D{E), which is roughly homogeneous, and the 
fractions of the eigenmodes a{E) in the initial condi- 
tion, which are almost equal. Hence, D{E) and a{E) are 
roughly constant[l3|. By setting u — exp(— /3i?)i, the in- 
tegral is rewritten as (1//3) (l/u)e~"du. By taking 
a limit of /? ^ oo first, \og{t) dependence is obtained 
asymptotically for large t. 

Though this estimate is originally asymptotic for large 
t, we used it for the time span where many eigenvalues 
contribute to the relaxation. For the last stage of the 
relaxation, only a few eigenvalues contribute. If there is 
a gap AA between two neighboring eigenvalues, there is 
a plateau in the relaxation for the time span For 
large (3, the gap between eigenvalues increases so that 
the existence of a plateau is expected. However, plateaus 
other than the last one, as well as their number during the 
relaxation, are not directly obtained from this argument. 
Here we give a heuristic argument for the plateaus. 

In Figs. 3 and 4, we give examples of the relaxation 
for smaller networks. Besides C{t), we have plotted 
"(t) = Xi{t)/xl'^ in Figs. 3(c)-(e) and 4(c). At each 
plateau, there are cluster(s) of elements in which xf™{t) 
takes almost the same value. Within each cluster, chem- 
icals are in local equilibrium through mutual reactions, 
whereas the equilibration process with elements out of 
the cluster is suppressed, since the concentrations of the 
catalytic components responsible for reactions for such 
equilibration are low. Consider a chemical with xf'^'" 
larger than the others. If the concentration of the cata- 
lyst (s) necessary to equilibrate the abundant chemical is 
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FIG. 3: (a) Time course of C{t) for two example networks, I 
and II, given in (b), with f3 = 16/e. In (b), the chemicals at- 
tached to the arrows that display reactions are their catalysts. 
Time courses of xf'^"{t) of (c) network I, and (d) network II 
corresponding to (a). In (a)(c)(d), all the chemicals are at 
equal number initially (i.e., {3 = 0). Time courses of xf'^'"{t) 
of network I from the initial condition Xi = 1 for i — 0,2, 3, 
xi = 0.4 and X4 = 1-6 are plotted in (e). 
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FIG. 4: (a) A reaction network with M = 12 and K = A, 
where the color of each arrow shows the catalyst for the re- 
action, and thick arrows indicate "major relaxation" for each 
chemical (see text), (b) Time course of C{t) for the network, 
and (c) time courses of xf^it) for /3 = 20/e. 



small, the equilibration process is suppressed. Negative 
correlation in the abundances between the excess chemi- 
cal and its catalyst will further suppress the relaxation to 
equilibrium. We now illustrate how this negative correla- 
tion gives rise to plateaus consisting of local-equilibrium 
clusters, by using examples given in Figs. 3 and 4. 

In the networks I and II in Fig. 3(b), consisting of 5 
chemicals, the component Xq (with lowest E) is trans- 
formed to all other components. In cases with large (3, 
because Eq is minimum, chemicals i > 1 flow into Xq 
from the initial condition with fi = 0, having xf'^'"{0) > 1 
for i > 1 for large /3. For both the networks, the eigenval- 
ues of J are exp(— /3i?i), exp{—(}E2), exp(— /3i?4), and 
{Ei — |e), asymptotically as /3 becomes large. As shown 
in Fig. 3(a), however, the numbers of plateaus appear- 
ing through the relaxation are different between the two 
networks. 

In network I, the first plateau consists of a local- 
equilibrium cluster Xq, X2, and Xi, whereas X3 joins 
to the cluster at the second plateau, as shown in Fig. 
3(c). The suppression of equilibration of Xi is explained 
as follows: Relaxation (i.e., decrease) of Xi (X4) is cat- 



alyzed by X4 (Xi), respectively. If one of the species 
Xi or X4 decreases faster, the relaxation of the other is 
suppressed. Because x^' is larger than x'^'' , X^ relaxes 
faster, so that the relaxation of Xi is suppressed. The 
negative correlation between the abundances of Xi and 
its catalyst hinders the relaxation of Xi . Since the relax- 
ations of X2 and Xi are catalyzed by the abundant Xi , 
the local-equilibrium Xq, X2, and X^ is first achieved 
and then X3 catalyzed by X2 (more abundant than Xi, 
the catalyst for Xi) joins to the cluster. 

In network II, on the other hand, the relaxation of Xi 
is not suppressed since its catalyst X2 relaxes only slowly 
because its catalyst X4 relaxes faster, as it is catalyzed by 
abundant Xi. Negative correlation exists, not between 
Xi and its catalyst, but instead between X2 and its cata- 
lyst Xi. Thus, the local equilibrium among Xq, Xi, X^, 
and Xi is realized to produce only one plateau. 

As expected from the above argument, the types of 
plateaus that appear in the relaxation can depend on 
the initial condition, because the reactions that are sup- 
pressed depend on which catalysts are first decreased. 
See Fig. 3(e), which shows the relaxation process of net- 
work I from the initial condition with Xi = 4x1. 
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For complex catalytic reaction networks, the argu- 
ment is not so simple, but the existence of local equi- 
libria and suppression of relaxation by the negative 
correlation mechanism generally underlie the formation 
of plateaus. Figure 4(a) is a catalytic reaction net- 
work with M = 12 and K ^ 4, and Fig. 4(b) (c) 
show the time courses of C{t) and xf^^{t) for (3 — 
20/e. As shown in Fig. 4(b), this network exhibits 
three plateaus in the relaxation process. At each 
plateau, chemicals i — 1, • • ■ M are clustered into a 
few groups within which xf™ is almost constant (Fig. 
4(c)). The following clusters are formed successively: 
{Xo,X4}, {X^jXq}, and {Xio, Xn} at the first plateau, 
{Xo,X3,X4,Xe}, {X^^Xr}, and {Xio,Xii} at the sec- 
ond plateau, and {Xq, Xi, • • • , Xy, Xg} and {Xio^Xu} 
at the third plateau. 

Each of these plateaus is explained by checking if the 
catalyst for the "majorh relaxation process for each Xi 
is abundant, which is a path catalyzed by Xk with the 
smallest k among the reactions with j smaller than i. At 
the first plateau, X^ and X4 have negative correlation 
since the major relaxation of X3 is the reaction catalyzed 
by X4, and that of X4 by ^"3, so that the formation of the 
cluster {Xq, X4} suppresses the equilibration between Xq 
and X3. 

At the second plateau, {X3, Xq} and {X5, X^} clusters 
have the following negative correlation. For {X^,X^} 
clusters, the reactions X^+Xq — > X2+XQ and X7-I-X3 
X2 -\- X'i give the major relaxations. On the other hand, 
the reactions Xq+Xj Xq+Xj and X3+X4 X-s+X4 
give the major relaxations for the {X3,Xq} cluster, but 
the reaction ^3 -I- ^"4 ^ X3 -I- X4 is suppressed since its 
catalyst X4 has already been decreased. In this case, the 
clusterjXs, X7} does not join X2, whereas the clusters 
{X3,Xe} and {Xo,X4} aggregate. 

In general, among a variety of chemical components, 
there exists such a negative correlation between chemi- 
cals in excess and the catalysts to decrease them towards 
equilibrium. Then, the equilibration of the chemicals is 
suppressed, leading to a plateau in the relaxation process. 

In this Letter, slow relaxation to equilibrium in cat- 
alytic reaction networks is demonstrated. When the tem- 
perature of the system is sufficiently lower than e/3, the 
average difference between energy levels, overall log(i) re- 
laxation appears. Several plateaus appear depending on 
the network and initial condition. The plateaus are not 
metastable states in the energy landscape but, rather, 
are a result of kinetic constraints due to a reaction bot- 
tleneck, originating in the formation of local-equilibrium 
clusters and suppression of equilibration by the negative 
correlation between an excess chemical and its catalyst. 

Possible configurations for local-equilibrium clusters 
are limited, and thus the number and ordering of plateaus 
are restricted. However, they are not necessarily uniquely 
determined by the network, but depend on the initial 
condition, because they are influenced by which cata- 



lysts are decreased first. Also, the relaxation is often 
non-monotonic; the deviation from equilibrium may in- 
crease during the relaxation course. Such roundabout 
relaxation has also been observed in a Hamiltonian sys- 
tem [III- We also note that discreteness in the molecule 
number results in anomalous reaction dynamics with long 
time correlations [31, and further sup pre sses the relax- 
ation in the catalytic reaction network |20l|. 

The behaviors reported here are reminiscent of the re- 
laxation in glass. Our model, as studied here, has a ki- 
netic constraint, although the constraint is based on the 
network structure rather than the spatial configuration. 
Application of theoretical frameworks developed in the 
study of glasses will be important to our chemical net 
glass in future work. Maintenance of the quasi-stationary 
states reported here, as well as successive changes in 
them, are often observed in biochemical processes, which 
have a large variance of reaction rates, i.e., potentiality 
of e » In future work, it will be important to dis- 
cuss the relevance of the present "glassyh dynamics to 
intracellular reactions. 
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